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Abstract 

The stability of nonlinear explicit difference schemes with not, in gen- 
eral, open domains of the scheme operators are studied. For the case of 
path-connected, bounded, and Lipschitz domains, we establish the notion 
that a multi-level nonlinear explicit scheme is stable iff (if and only if) 
the corresponding scheme in variations is stable. A new modification of 
the central Lax-Friedrichs (LxF) scheme is developed to be of the second 
order accuracy. The modified scheme is based on nonstaggered grids. A 
monotone piecewise cubic interpolation is used in the central scheme to 
give an accurate approximation for the model in question. The stability of 
the modified scheme is investigated. Some versions of the modified scheme 
are tested on several conservation laws, and the scheme is found to be ac- 
curate and robust. As applied to hyperbolic conservation laws with, in 
general, stiff source terms, it is constructed a second order nonstaggered 
central scheme based on operator-splitting techniques. 

1 Introduction 

We are mainly concerned with the stability of nonlinear explicit difference 
schemes arising, e.g., in numerical analysis of nonlinear PDE (partial differ- 
ential equation, such as parabolic, hyperbolic, etc.) systems, generally, in mul- 
tidimensional space. Stability is the central and the most pressing problem in 
any algorithm [35]. Nowadays, there exists a few methods for stability analysis 
of some classes of difference schemes approximating PDE systems (see, e.g., [16], 
[17], [32], [35], [38], [48] and references therein). However, the problem of sta- 
bility analysis for nonlinear schemes is still one of the most burning problems. 
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because of the absence of its complete solution [T^- In particular, the vast ma- 
jority of difference schemes, currently in use, have still not been analyzed [16], 
and, in general, no numerical method for non-linear systems of equations has 
been proven to be stable [HI]- A new approach to testing scheme stability has 
been suggested in [7] to prove convergence of non-linear schemes for systems 
of PDEs. It was demonstrated [7] that the notion of scheme in variations (or 
variational scheme [9], [10]) has much potential to be an effective tool for study- 
ing stability of nonlinear schemes. The scheme in variations for a difference 
scheme represents a tangent space at a point of the manifold associated to the 
difference scheme [7], [S]- Thus, the scheme in variations will always be linear 
and, hence, enables the investigation of the stability for nonlinear operators us- 
ing, in general, an infinite family of linear patterns. It was established in [7] 
the notion that the stability of a scheme in variations implies the stability of 
its original scheme, and that a nonlinear explicit scheme will be stable ijf (if 
and only if) its scheme in variations will be stable. In this paper we consider 
difference schemes with not, in general, open domains of the scheme operators. 
For the case of internal path-connected, bounded, and Lipschitz domains, we 
establish the notion that a multi-level nonlinear explicit scheme is stable iff the 
corresponding scheme in variations is stable. 

By way of illustration of the developed approach we are concerned with the 
stability analysis of central difference schemes for hyperbolic systems of conser- 
vation laws with source terms. Such systems are used to describe many physical 
problems of great practical importance in magneto-hydrodynamics, kinetic the- 
ory of rarefied gases, linear and nonlinear waves, viscoelasticity, multi-phase 
flows and phase transitions, shallow waters, etc. (see, e.g., [6], [12], [18], [24], 
p9] . [32] . [34] . [37] . [41], [42]). We will consider a system of hyperbohc conser- 
vation laws written as follows (e.g., [TH], [32) 

^+E5^f'(") = 7q(")' 0<i<?lnax, u(x,t)|,^o = u"(x), (1) 

where the column- vector x = {xi,X2, ■ ■ ■ ,xn}'^ G M^, u — {ui, U2, ■ ■ ■ , mm}"^ is 
a vector- valued function from x [0, -foo) into a subset fiu C M*^, ij (u) = 
{fij (u) , f2j (u) , . . . , fMj (u)}"^ is a smooth function (flux- function) from fiu 
into M.^ , q (u) = {qi (u) , (72 (u) , . . . ,qM (u)} denotes the source term, r > 
denotes the stiffness parameter, u'^ (x) is of compact support. It is assumed 
that T = const without loss of generality. It is also assumed that all eigenvalues, 
— (u) , of the Jacobian matrix G (=9q(u) /9u) have non-positive real 
parts, i.e. 

Re^k (u) < 0, Vfc, Vu ef^u. (2) 
In what follows |lM|jp denotes the matrix norm of a matrix M induced by the 

vector norm ||v||p = (J^i biD^^^'i and ||M|| denotes the matrix norm induced 
by a prescribed vector norm. M denotes the field of real numbers. The transpose 
of a matrix (or vector) S is denoted by S^. The null element in any linear space, 
as well as the number zero, will be denoted by the same symbol 0. 
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Central schemes are attractive for various reasons: no Riemann solvers, char- 
acteristic decompositions, complicated flux splittings, etc., must be involved in 
construction of a central scheme (see, e.g., [S], [5], (33], [35], [3^, [HI], [31], 01], 
[33] and references therein). However, there is the need for staggered central 
schemes to alternate between two staggered grids, which can be cumbersome, 
e.g., near the boundaries [33]. Moreover, there is a risk for every central scheme 
to exhibit spurious solutions [8] in spite of sufficiently small CFL (Courant- 
Friedrichs-Lewy |32j ) number. To simplify implementation of central schemes 
in the case of complex geometries and boundary conditions, it was developed 
nonstaggered central schemes in [33]. A new second-order modification of the 
staggered LxF scheme was developed in [5] to avoid the risk of spurious oscilla- 
tions. It was demonstrated in [8] that the higher order versions of LxF scheme 
can produce spurious oscillations because of a negative numerical viscosity in- 
troduced to increase accuracy of the scheme up to 0(At -|- (Ax)^). To reduce 
the risk of spurious solutions, it was suggested in [S] to introduce an additional 
non-negative numerical viscosity such that the scheme's order of accuracy be- 
comes 0{{At)^ + (Aa;)^). The developed central scheme was tested on several 
conservation laws taking a CFL number equal or close to unity, and the scheme 
was found to be accurate and robust. In this paper, we extend the second-order 
modification to a new nonstaggered central scheme. The stability of this scheme 
is proven in Section 13.31 on the basis of the approach developed in Section [3] as 
well as in [7], [5]. The scheme is tested on several conservation laws in Section 

El 

A stable numerical scheme may yield spurious results when applied to a 
stiff hyperbolic system with relaxation (see, e.g., [1], [4], [6], [12], [13], [24] . 
[44] . [45]). It is significant that a numerical scheme for relaxation systems 
must possess a discrete analogy to the continuous asymptotic limit, because any 
scheme violating the correct asymptotic limit leads to spurious or poor solutions 
(see, e.g., [12], [24], [25], [37], [41]). Most methods for solving such systems can 
be described as operator splitting ones, [T3], or methods of fractional steps, 
[B] . After operator splitting, one solves the advection homogeneous system, and 
then the ordinary differential equations associated with the source terms. We 
are mainly concerned with such an approach in Section UJ 



2 Stability of nonlinear explicit schemes 

We consider the following {q + l)-level difference scheme: 

w7,w7-i,...,wr'+'), :r!„CM^^M^°, (3) 

where q, Nq, I > 1 arc finite integer constants, N — qNoI, E denotes a 
vector-valued grid function, i G uJi denotes a node of the grid cui = {1, 2, . . . , /}, 
k E UJ2 denotes a node (time level) of the grid uj2 = {0, 1, ... , M}, M > q is a 
finite integer constant. G"= {C,"!, G-'jj ■ • • , No) denotes the vector- valued 
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function with the domain and range belonging to M.^ and , respectively. We 
emphasize that the domain, r2„, of the function Gf is the set of argument values, 
which need not be open in IR^. We assume that n in (jS]) denotes the time level, 
tn (= nAt). Thus, the time increment will be represented by At = tma^/M = 
const, where ^max denotes some finite time over which we wish to compute. 
Using the notation 




we reduce (see, e.g., [IS], [37]) the multi-level scheme, to the following 
two-level scheme: 

v,"+i =Hnvr,vJ,...,v7), :1]„CM^^M«^«, lewi, n,n+lew2. (5) 
If we define 

v"= {K)^ , K)^ , . . . , (v?)^}^ , H"= {(H^)"^ , (H«)^ , . . . , {WifY , (6) 

then the scheme (O becomes 

v"+i=H"(v"), H":aiCR^^R^, n, n + 1 e = {0, 1, . . . , M} . (7) 

As usual (e.g., gOl p. 62]), for mappings f : 17/ C ^ and g : 17g C 
— >■ R^, the composite mapping h = g o f is defined by h (v) = g (f (v)) for 

all V sri/i = {v gf]/ I f (v) € ilg}. Using the composite mapping approach, we 

rewrite Scheme ([7]) to read 

y = F (x) , F : 17f C R^ ^ R^, (8) 

where the following notation is used: x = v", y = v^^, F = H*^~'^oH*^^^o. . .o 
HO, = {v"el}o I vi = HO (vO) e r!i I . . . I v*^-i = H*^-^ (v^'^-^) g ^m-i). 
Notice, F in ([8]) depends also on scheme parameters (e.g., space and time incre- 
ments), however, this dependence is usually not included in the notation. Let 
the scheme parameters (including time increments) be represented by a vector 
s belonging to some normed space with the norm |s|. 

Let us remind the well-known definition of stability of an explicit scheme 
(see, e.g., [16], [18], [32], [47], [48], [49]). 

Definition 1 Scheme ^Bj) is said to be stable if there exist positive constants sq 
and C such that for all x, x, € fi^? the following inequality is valid 

|lF(x,)-F(x)|| <C||x, -x||, Vs:|s|<so. (9) 
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Thus, scheme ([5]) and, hence, scheme ^ wih be stable ijf the function F will 
be Lipschitz for a constant C. Notice that this well-known definition of stability 
may lead to wrong conclusion. Actually, let us consider the "slit plane" [3T] in 
polar coordinates (r, 9) 

np = {(r, e)\0<r <oo, -TT < < tt} C (10) 

and the vector- valued function F — {Fi (r, 9) , F2 (r, 9)}'^ such that [51] 

Fi = r, F2 = 9/2. (11) 

We take (r, 9)^ = (rg, — tt -|- e), (r, 6*) = (ro, tt — e), and ro = const. Considering 
the Euclidean distance between these two points of the plane with Cartesian 
coordinates = (—tq cose, —rg sine) and x = (—tq cos e, rg sine) we find that 
the inequality © takes the form: cos(e/2) < Csine. Obviously the mapping 
PTj) is not Lipschitz, since C in ^ tends to infinity as e — > 0. Therefore, we have 
to conclude, in view of the above definition, that the scheme ([TT|) is not stable, 
even though for every point {r,9) in Hp, (ITU)) , there exists a neighborhood of 
(r, 9) such that the function F, pT|) . restricted to the neighborhood is Lipschitz 
for C — I, i.e. the function F is locally Lipschitz. Further, F is the non- 
stretching mapping of the "slit plane" pUj) into the right semi-plane. Hence, the 
preceding definition of stability, i.e. Definition [U needs to be slightly improved. 
Such a problem has been considered in [7] under the assumption that the domain 

of the vector- valued function F, ([5]), is open in R^. We will assume that 
the Lebesgue- measurable domain fi^ of the function F, ([81), is a not, in general, 
open subset of M.^ with non-empty interior. 

We shall start with some basic notions. A set fl C R^ is said to be internal 
path-connected (or path-connected, in the case that ft is open) if every two 
points X, X* G O can be joined by a continuous curve 7 : [0, 1] C R — of finite 
length, L{j), with int{'j) C int{n). Here and in what follows int {fl) denotes 
the interior of C R^ and int (7) = 7\(7 (0) U 7 (1)) denotes the interior of 
the continuous curve 7. The intrinsic metric on an internal path-connected 
set fl is defined as 

An (x, X,) = inf L (7) , x = 7 (0) , x, = 7 (1) , x, x,el^. (12) 

7Cr2 

An open ball (of radius r > 0) about x eR^ is denoted by B (x, r) (or just B-x). 
A closed ball is denoted by B (x, r) (or just B^). Given points x, y SR^, given 
an open ball _B(x, Tx), and given an open ball B{y,ry) such that x ^By, the 
convex open set Qx = {z | x + A (z — x) , A > 0, z GBy} n B^ will be called 
a finite cone with vertex x. A set O C R^ is said to have the cone property 
if there exists a finite cone such that each boundary point x € dH, is the 
vertex of a finite cone Qx ^ int{fl) congruent to Q^. Suppose that il C R^ 
has the cone property, and let denote the union of all finite cones with the 
vertex x edQ, i.e. = UQx- Then a function F : C R^ -> R^ wiU be 
called locally Lipschitz (or locally Lipschitz continuous) if for every x Gfi there 
exist _Bx and a constant Cx > such that \\F (x) — F (y)|| < Cx ||x — y|| for each 
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y G-Bx provided x Eint (fl), or y GBxI^K^ otherwise, ft C is Lipschitz (has 
a Lipschitz boundary), if for each point x GdU there exists an open ball such 
that dil n Bx is the graph of a Lipschitz function and int (fi) lies on one side of 
its boundary (for more details see, e.g., [TTJ p. 149], [31] p. 354]). Let us note 
that ri C satisfies the cone property if it has a Lipschitz boundary (see, 
e-g-, [HI P- 151], [31, p. 355]), but not vice-versa as the example (ITU)) shows. 
W^'°° (fi) denotes the Sobolev space (see Definition 5 in [TTJ p. 28], see also [21] 
and [3T]), i.e., the space consisting of all bounded functions on open fl C 
whose distributional derivatives are bounded functions as well. Let F = {Fi , 
F2, Fn}^ in ([U, let VFi = {diFi, dNFi} denote the distributional 
gradient of Fi, and let 5F, Sx € denote variations. The following equality 

(5F = F'.5x, F'= {(V^^i)^,(VF2)^,...,(Vi^^Ar)^}^, (13) 

will be viewed as the scheme in variations for ([5]). The matrix representation 
of F' in ([T3)) is given by the following Jacobian matrix [TD]: F' = {dFi/dxj}, 
i,j ^1,2,... ,N. 

Definition 2 Let flp in be internal path- connected. Scheme is said 
to be stable if there exist positive constants sq and C such that the following 
inequality holds 

||F(x,)-F(x)|| <CAo^(x,x,), Vx,x,g17f, V s : |s| < sq- (14) 

Notice that the scheme PT|) is stable in terms of Definition [51 since the 
inequality (|14p holds for C = 1. It is significant that the domain Vtp, (|10l) . of 
the function F, (|lip . is open in R^. However, if the domain of the function F, 
pip , is not open, then the situation can be in complete contrast to the previous 
one. Actually, let us add to Vtp, (ITO)) . the only boundary point (l,7r), i.e. let 

= r2i^ U (1, tt) will be the domain of the function F* = F, ([TT]), if (r, 6) eVIf 
andF* = (l,7r/2) if (r, 6*) = (1, tt). We take (r,6l), = (l,-7r + e), (r,^) = (l,7r). 
Considering the Euclidean distance between these two points of the plane with 
Cartesian coordinates x* = (— cose, — ro sine) and x = (—1, 0) we find that the 
inequality p4|) takes the form: cos (e/4) < C sin (e/2). Thus, the scheme (ITTj) . 
namely F : Wp C — >■ M^, is not stable, since C in (|T4)) tends to infinity as 
e ^- 0. 

Lemma 3 Let the path- connected Hp of (0) be open in M^. Scheme (0) will be 
stable in terms of Definition\^iS F in ^ will be locally Lipschitz for a common 
constant C , for all scheme parameters s such that |s| < sq. 

Proof. Suppose Scheme ([8]) is stable, i.e. is valid. Choose any point 
xef^F- Since fi^? is open, there exists a radius r such that B{x.,r) C ftp. 
Choose any point x, e B (x, r) , and let 7^ be the straight line segment joining 
the points x, x* € i? (x, r). In view of p^ . F in ([5]) will be locally Lipschitz for a 
common constant C, for all s : |s| < sq, since A^^ (x, x,) = L (7^) = ||x, — x||. 
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Conversely, suppose that F in ([5|) is locally Lipschitz for a common constant 
C, for all s : |s| < Sq. Let some points x, x* g be joined by a continuous 
curve 7. In view of ([T^ , the curve 7 can be taken such that L (7) < Aq^ (x, x* )+ 
e for an arbitrary e > 0. Given any point z S 7, there is a ball B^, C f2_F 
such that F restricted to is Lipschitz for the common constant C (locally 
Lipschitz continuity). The balls {B^} form an open cover of 7. Since the 
mapping 7 : [0, 1] C M — )■ is continuous, the curve 7 is compact [28l p. 94]. 
Hence, by the compactness of 7, {B^} has a finite subcover consisting of balls 
-Bx = Szi, ^22, . . ., B^fi = ^x*- Since F is locally Lipschitz, we find 

||F(zfc+i) -F(zfc)|| < C7||zfc+i -Zfcll , fc = l,2,...i^-l, Vs: |s| < sq. (15) 

Then, by virtue of p^ . we find 



|F(x.)-F(x)| 



^[F (zfc+i)-F(zfc) 



V 

k 

CL{^) <C An A^,^*)+eC, V s : |s| < sq- (16) 



< C ^ ||zfc+i - Zfcll < 

k 



By letting e — > 0, we find that ([T4|) holds. ■ 

Notice, in view of Lemma[3l the scheme (fTT]) is stable in terms of Definition[2l 
since the domain ftp, (|T0| . is open in and path-connected, and the function 
F, dTTI) . is locally Lipschitz for the common constant C — 1. 

Theorem 4 Let the internal path- connected ftp of ^ have the cone property, 
then the scheme (0) will be stable in terms of Definition iff F in ^ will he 
locally Lipschitz for a common constant C , for all s such that |s| < sq. 

Proof. Let denote the union of all finite cones with the vertex x ^dflp. 
Suppose Scheme ([5]) is stable. Choose any point x gJIf- If x; £int (ftp), then, 
in view of Lemma [3l F in ([S]) will be locally Lipschitz. If x ^int (Hh), then 
(cone property) for each x» G there is a straight line segment joining the 
points X, x». Since Aq^ (x, x») — ||x, — x||, the necessity, in view of p^ . is 
proven. 

Conversely, suppose F in is locally Lipschitz. Let some points x, x, S fi^ 
be joined by a continuous curve 7. If x, x, e int (^p), then, in view of Lemma 
[3l we find that (fT4|) holds. Let s. ^int (yip), however x, e int (ftp). Choose 
any point z G Kx such that ||x — z|| = e. Then, since F is locally Lipschitz, we 
have IjF (x) — F (z)|| < C ||x — z|| for a sufficiently small e. Since z e int (Hp), 
we find, in view of Lemma [31 that ||F (z) — F (x*)|| < CAfip (z,x»). Then 

||F (x) - F (xOII < ||F (x) - F (z)|| + ||F (z) - F (x,)|| 

<eC+||F(z)-F(x*)||, Vz: ||x-z|| =£, Vs: |s|<so- (17) 
Since ||F (x) — F (x,)|| does not depend on z, we find 

||F(x)-F(xO|| < inf {eC+||F(z)-F(x,)||} (18) 

z: ||x— z||— e 
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<eC + C inf An^(z,x*), Vs: |s| < sq. (19) 

z: ||x— z||— e 

By letting e — ?> 0, i.e. z — >■ x, we find that (fH)) holds. The proof in the case x, 
X* ^ int {^f) is reduced to the previous one by choosing z g i^x, • ■ 

We shall need the following theorem that identifies the Sobolev space 14^^'°° [Vl) 
as a space of locally Lipschitz functions. 

Theorem 5 Let ^Ip C he open, and let F = {Fi , F2, . . . , Fn}^ in 
Then, Fi, i = 1,2,...,A^, (and, hence, F) is locally Lipschitz (in the sense of 
having representatives) iff Fi G W^'°° {^f)- 

The proof of Theorem [5] can be found, e.g., in Theorem 4.1], [211 p. 
342]. The following lemma gives the necessary and sufficient conditions for the 
stability of the scheme in variations (I13p . 

Lemma 6 Linear Scheme US\) will be stable iff there exist positive sq,C — const 
such that 

\\F'\\<C = const, Vxef^F, Vs: |s| < sq. (20) 

The proof of Lcmma[6]can be found in (see also, e.g., [16], [47], [48], [49]). 
Let us remind the well-known theorem that ensures that Lipschitz domains are 
extension domains (see ^31, pp. 320, 356], [iTJ pp. 149, 285]). 

Theorem 7 Let Vt d be an open set with uniformly Lipschitz boundary. 
Then there exists a continuous linear operator £ .■ W^'°° (f2) — > (M^) 
such that £(u)(x) = m(x) for all u S W'^^°^ {n). 

The interconnection between stability of the nonlinear scheme (jSj) and sta- 
bility of its scheme in variations is established by the following theorem. It is 
assumed that fi^ in ([8]) need not be open. It is also assumed that the set fii? is 
bounded, i.e. there exists rg < 00 such that C Bq = B (0, rg). Let F = {Fi , 
F2, . . . , Fn} and let VF^, i = 1, 2, . . . , iV, denote the distributional gradient 
of Fi. Definition of the distributional gradient (or distributional partial deriva- 
tive) of at X G int{n,p) can be found, e.g., in [3T]. If x £ d^lp, while Hp 
is Lipschitz, then the distributional gradient, VF^, is assumed to be equal to 
V£(F,)(x). 

Theorem 8 Consider the scheme (0). Let F be bounded and let Qp be internal 
path- connected, bounded, and Lipschitz. Then, the scheme, (0), will he stable in 
terms of Definition\^iS its scheme in variations, I113\) . will be stable. 

Proof. Let the scheme ^ be stable, then, in view of Theorem HJ F will 
be locally Lipschitz for a common constant C, and, hence, in view of Theorem 
[3 F, e {int {Qp)), i ^ 1,2,...,N. By virtue of Theorem [7] we can find 

F* e (M^) such that i^*(x) = Fi(x) onintiflp). Hence F* e W^'°° {Bq), 

where Bq is an open ball such that Hp c Bq c R^. Notice, ||Vi^*||^, i = 
1,2, . . . , N , is bounded on the ball Bq and, naturally, on the domain flp d Bq, 
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since F* G W^ °° {Bq). In such a case there exists a common constant C such 
that IIF'II < C for all x g Qp, where ||F'|| = ||f;^||, = {||VFi||^ , IIVF2IL, 
. . . , II Vi^jv I loci Hence, in view of Lemma [BJ the scheme in variations, (IT^ . is 
stable. 

Conversely, suppose that the scheme in variations, (|13l) . is stable. In view 
of Lemma m there exists a common constant C such that ||F'|| < C on flp, 
and, hence, HVFJ^, i = 1,2,...,A^, is bounded on V,p. Consequently, Fi S 
W^'°° {int {^Ip)), i = 1,2, ...,A^. Since intiVlp) is an extension domain for 
W^'°° {int {np)), we can find, in view of Theorem H F* e W^-"^ (M^) such 
that i^*(x) = i^i(x) on int{np). In view of Theorem[5l F* , i ^ \,2, . . . ,N , wih 
be locally Lipschitz on (and, hence, on Vlp) for a common constant. Then 
Fi, i = 1,2, . . . , N, and hence F will be locally Lipschitz on il^, since Hp being 
Lipschitz satisfies the cone property. Thus, in view of Theorem |4l the scheme 
(O wiU be stable. ■ 

Notice, if F in the scheme ^ is Gateaux-differentiable, then VF^ (see The- 
orem |8]) denotes the classical gradient of Fi, and, hence, it may be taken that 
fi. = {||VFi|| , IIVF2II, IIVF^Ilf . 

Remark 9 Definitions^ is "restored in its rights" for domains being internal 
path- connected, bounded, and Lipschitz (see the proof of Theorem\^. Namely, 
considering the scheme with such a domain, we find that the scheme will be 
stable in terms of Definition\^iS it will be stable in terms of Definition\^ 



3 Construction of stable central schemes 

In this section we consider explicit schemes on a uniform grid with time step Ai 
and spatial mesh size Ax, as applied to 1-D hyperbolic equations. In view of 
the CFL condition (32], we assume for the explicit schemes that At — O (Ax). 
Moreover, we will also assume that Ax = O (At), since a central scheme gen- 
erates, in general, a conditional approximation (see [8]). In such a case, the 
following inequalities will be valid, for sufficiently small At and Ax, 

voAt < Ax < ^if^At, i^OjMo ~ const, Q < < /Xq. (21) 

We will focus on the following 1-D version of the problem ([T]): 

^ + — f (U) =0, tn<t< tn+l = tn + At, U {X, t„) = U" (x) , (22) 

Since the system ((22|) is hyperbolic, the Jacobian matrix of f (u) possesses M 
linearly independent eigenvectors (see, e.g., [Hj). In addition, it is also assumed 
that 

sup IIAII2 < Aniax < 00, A =^^. (23) 

Using the central differencing, we write 



t = t„ + 0.25, X = Xi^.o.5 



dt 
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£n+0.25 _ fn+0.25 



dx 



^+1 



Ax 



o 



(i^xf) . (25) 



By virtue of we approximate (1^^ on the cell [xi,Xi+i\ x [i„,t„+o.5] by 

the following difference equation 

n+0.5 _ n _ ( n+0.25 _ n+0.25\ /2C^ 

As usual, the mathematical treatment for the second step (i.e., on the cell 
[a;i_o.5, Xi+o.s] X [t„+o.5, tn+i]) of a staggered scheme will, in general, not be 
included in the text, because it is quite similar to the treatment for the first 
step. 

Considering that approximates ([^ with the accuracy 0((Aa;)'^+(At)^), 
the next problem is to approximate vJYq g and g"^*^'^^ in such a way as to retain 
the accuracy of the approximation. For instance, the following approximations 

vr+o.5=0.5(v»+vr+0+o((Ax)2), g,"+"-25=f(vn+0(Ai), (27) 

leads to the staggered form of the famed LxF scheme that is of the first-order 
approximation (see, e.g., [HI p. 170]). One way to obtain a higher-order scheme 
is to use a higher order interpolation. At the same time it is required of the 
interpolant to be monotonicity preserving. Notice, the classic cubic spline |46) 
does not possess such a property (see Figure [T^). In the following, subsection 
13.11 we consider the problem of high-order interpolation of v"^q g in (|26p with 
closer inspection. 

3.1 Monotone piecewise cubics in construction of ex- 
plicit schemes 

We will consider some theoretical aspects for high-order interpolation and em- 
ployment of monotone piecewise cubics (see, e.g., [14], [27]) in construction 
of monotone explicit schemes. 

Let p = p {x) = {p^ {x) , . . . ,p'' (x) , . . . (a;)} be an interpolant, and let 

Pi = P (a^j) , P't ^ P' (xi) , Ap, = Pj+1 - p,, 

P^=A..^, P'...=B..^, (28) 

where p^ denotes the derivative of the interpolant at x = Xi. The diagonal 
matrices and in (|28l) are defined as follows 



A, = dtag {alal . . . , a^} , M, = diag ...,PT]- (29) 

If p = p {x) is a piecewise cubic interpolant (see, e.g., [T3], [27]), then it will 
be component-wise monotone on [xi,a;i-|-i] iff one of the following conditions 
(see [H], [27]) is satisfied: 

{al - if + (af - 1) - l) + - l)' - 3 (al + /?f - 2) < 0, (30) 



10 



3, a?>0, /3?>0, Vi,fc. (31) 

The region of monotonicity is shown in Figure [T}). The resuhs of implementing 
a monotone piecewise cubic interpolation when compared with the classic 
cubic spline interpolation, are depicted in Figure [T^. 




0123456789 10 1 2 3 4 



Figure 1: Monotone piecewise cubic interpolation, (a) Interpolation of a 1-D 
tabulated function. Circles: prescribed tabulated values; Dashed line: classic 
cubic spline; Solid line: monotone piecewise cubic under a, /3 < 3; Dotted line: 
monotone piecewise cubic under a,/3 < 1. (b) Necessary and sufficient condi- 
tions for monotonicity. Horizontal hatching: region of monotonicity; Unshaded: 
cubic is non-monotone. 

We note (Figure [T^) that the constructed function produces monotone in- 
terpolation and has a small, or even practically zero, deviation from the classic 
cubic spline at some sections where the classic cubic spline is monotone. 

Let us consider the problem of monotone high-order approximation of v"_|_q g 
in (f26|) . First, given the values, Pi and Pi+i, of the interpolant p — p{x) and 
the approximate estimates, and d^+i, of its derivatives, namely 

p^ = d, + 0((Ax)^), (32) 
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we obtain the following interpolation formula 



p,+o.5 = 0.5 (p, + p,+i) - — (d,+i - d,) + O {{AxY) , (33) 

where 

r = min(4, s + 1) . (34) 

Actually, let p (x) be sufficiently smooth, then, using the Taylor series expansion 
of the function p (x), we obtain: 

p,+i + p, = 2p,+o.5 + p:'+05 (^)' + ^ (^)' + ^ ((^^)') ■ (^^^ 

In a similar manner, using the Taylor series expansion of the function p' (x), we 
can show that 

P^+i - p: = p'/+05Ax + ^ [^j + O [{Axf) (36) 
By virtue of (1551) we obtain from (1551) that 

p,+o.5 = 0.5 (p. + p,+i) - ^ (p',+1 - pO + ^ (^) ' + O [{Ax)') (37) 
In view of p2p we obtain from ([37]) the following interpolation formula 

p,+o.5 = 0.5 (p, + p,+i) - ^ (d,+i - d,) + O [{Axf + {AxY+^) . (38) 

Hence, by virtue of ([38]) and ([33]) . we conclude that ([34| is, in general, valid. It 
is clear that if ([32]) is valid and d^+i — d^ = (Pi+i — Pi) + 0{{AxY'^^), then 
r = min (4, s + 2) in ([55]) . Thus, if p (x) has a continuous fourth derivative and 
p^ can be estimated with accuracy OiiAxf) (see, e.g., [H p. 112]), then r = 4 
in ([33|. 

Let us note that instead of point values employed in the construction of the 
scheme ([26]) . it can be used the cell averages (see, e.g., [5], [29], |32]). Then, 
given the values, Pj, of the averaged interpolant p = p(a;) and the estimates, 
di and d^+i, of its derivatives, namely p ■ = d^ + O {{AxY) we obtain the same 
formula ([551) written for the averaged values. Further, to estimate the flux 
at Xi, the point value Pi can be reconstructed using a piecewise polynomial 
interpolation (e.g., |S], [22] )■ Notice, we need not any reconstruction when 
dealing with first- and second-order schemes, since the point values agree with 
the corresponding cell averages to 0((Ax)^) (e.g., [29], [32]). We will therefore 
omit the bar notation when dealing with such schemes. 

Instead of ([55]) . we will employ the following interpolation formula 

Ax 

Pi+0.5 = 0.5 (pi -f Pi+i) - ><— (di+i - di) , (39) 
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where, for some reason, a new parameter >c is used. Consider, for instance, the 
case when the point value vf^_Q 5 in the scheme ([26| is estimated by the cell 
average calculated on the basis of the monotone piecewise cubics. Note that 
this estimation is of second-order. Such an approach leads to the interpolation 
of the point value by formula (15^]) with the parameter yt = 2/3. Obviously, 
the interpolation ([55)1 is of second-order under < x < 1, and ([55)1 coincides 
with p3|) under x — 1 to give fourth order accuracy. Notice, if < < 1, 
then Pi+0.5 can be estimated with accuracy 0((Aa;)^) by provided that 
1 — X = 0{{Ax) ). Actually, by virtue of ([55)) and (15^ . we obtain that 

K— (d,+i - dO = — (d.+i - d,)+{l^>c y '+^^ = — (d,+i - dO 

+ + o {{Axr') ^ ^ (d.+i - d,) + o{{Axn m 

The approximation of derivatives can be done by the following three 
steps [13]: (i) initialization of the derivatives = d^ -f O {{AxY); (ii) choice 
of subregion of monotonicity; (iii) modification of the initialized derivatives, d^, 
to produce a monotone interpolant. 

The matter of initialization of the derivatives is the most subtle issue of 
this algorithm. Thus, using the two-point or the three-point (centered) differ- 
ence formula (e.g. [27], [41]) we obtain, in general, second-order approximation. 
Performing the initialization of the derivatives p^ = d^ -I- O ((Ax)*) in the in- 
terpolation formula ([39)) by the classic cubic spline [46] interpolation, we obtain 
s = 3 (e.g., [1^, [17]). The same accuracy can be achieved by using the four- 
point approximation [271. However, the efficiency of the algorithm based on 
the classic cubic spline interpolation is comparable with the one based on the 
four-point approximation, as the number of multiplications and divisions per 
one node is approximately the same for both algorithms. 

Obviously, for each interval [xi,Xi^i] in which the initialized derivatives d^, 
di+i such that at least one point (af, /3*^) does not belong to the region of 
monotonicity ([50 )) - ([5T)) . the derivatives d^, d^+i must be modified to d^, d^+i 
such that the point (5^ , /3,j ) will be in the region of monotonicity. Moreover, it 
is desirable to modify the derivatives d = {. . . , d^;^, df , d^^, ■ ■ in such a 

way that d = |. . . , d^^, df , d^^, . . .| would be the solution to the following 
mathematical programming problem: 



d-d 



niin. (41) 

d 



The modification of the initialized derivatives, would be much simplified if we 
take a square as a subregion of monotonicity. In connection with this, we will 
make use the subregions of monotonicity represented in the following form: 

< a,^ < 4H, 0<P-< 4H, Vi, k, (42) 
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where H is a monotonicity parameter. If H = 0.75, then the subregion (H^ , called 
de Boor-Swartz box [27] , coincides with the one used by Fritsch and Carlson [14] . 
It is proven in [7] that the interpolation (P^j) will be monotone, i.e. the value of 
an arbitrary component of Pi+o.5 will be between the corresponding components 
of Pi and Pi+i, iff dHJ) will be valid provided that < H < 1. Notice, in such 
a case the point (aj^, /3^) may be out of the region of monotonicity shown in 
Figure [1)3. 

To fulfill the conditions of monotonicity the modification of derivatives 
di = [djjdf, . . . , d™} can be done by the following algorithm suggested, in fact, 
by Fritsch and Carlson [TJ] (see also [17]): 



:= 4H minmod(Af_i, A,^), := minmod(df , ), H = const, (43) 



3.2 Construction of nonstaggered central schemes 

We will consider explicit central schemes on a uniform grid with time step At and 
spatial mesh size Ax. In view of the interpolation formula (1391) . the staggered 
scheme (|26|) becomes 

- 0.5 (v^Vi + vD - (d^Vi - dD - T^^^H^^' (4^) 

where d" denotes the derivative of the interpolant at x = Xi, the range of values 
for the parameter x is the segment {) < k < 1. li x — 0, then the scheme 
(|45|) coincides with the LxF scheme. It is shown in [^ that the first-order, 
0{At+{AxY /'At+{Ax) ), scheme (|^5|) generates a conditional approximation, 
because it approximates (1^^ only if {AxY /Ai — )■ as Ace — )■ and At — )■ 0, 
where r is the order of approximation of v"^g 5 by the interpolation formula 
(jgg)) . The scheme ^ is abbreviated in [5] as COSl. 

Let us develop new nonoscillatory central schemes, which are based on regu- 
lar, nonstaggered spatial grids. Using the central differencing, we approximate 
(|22|) on the cell [xi^i, Xi+i] x [tn,tn+i\ by the following difference equation 



Notice, the approximation g"^°'^ — f (v") -I- 0{At), leads to the first-order, 
0{At+{Axf), scheme 





(46) 




(47) 
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which is absolutely unstable [111 P- 113], [321 P- 71]. Let us construct a stable 
scheme on the basis of (|46)) . The interpolation formula (p9)) . as applied to the 
cell [xi-i,Xi+i\ X [tn,tn+i\, becomcs 

Aa; 

= 0.5 (p.j_i + pi+i) - >f— (di+i - di-i) . (48) 



Equality ([48]) . after elementary transformations, leads to 

Ax 

p, = 0.25 (p,_i + 2p, + p,+i) - X— (d,+i - di_i) . (49) 



Using (|49l) we obtain from (l47l) the following nonstaggered central scheme 

where < >f < 1. In view of as applied to the cell [a;i_i,Xi+i] x 

[tm ^n+i], the local truncation error, on a sufficiently smooth solution u (x, t) 
to (HH) is found to be 



2 ' 



^ = 0(At + (Ax)^) + O {^-^^ + (1 - x)0 (^i^ j , (51) 

where r = min(4, s + 1), s denotes the order of approximation in ([32]). Note 
that the scheme (|50| will be 0{At + (Ax) ) accurate if r > 3 as well as 1 — >f = 
0(Ax). The first-order nonstaggered central scheme ([50]) . being a monotone 
approximation of the 1-D equation ([22]) by cubics, will be abbreviated to 
MACl. 

Let us consider the transformation from the scheme (|Tf| to ([501 with closer 
inspection. We will use the, so called, first differential approximation of the 
scheme ([47)) ([151 P- 45], [49l p. 376]; see also 'modified equations' in [151 P- 
45], [21], [Sij). As reported in [TS], [H], this heuristic method was originally 
presented by Hirt (1968) (see [TSl p. 45]) as well as by Shokin and Yanenko 
(1968) (see p. 376]), and has since been widely employed in the development 
of stable difference schemes for PDEs. 

The first differential approximation of the scheme ([47]) is the following: 

The negative diffusion term in ([5^ is associated with instability, as this therm 
is the source of energy resulting in an unlimited growth of the amplitude of the 
solution (see, e.g., [FT], [S5]). Notice, for the sake of convenience we use the 
same notation in ([^^ and in ([5^ in spite of the fact that these equations are 
different. We will use such an approach if it does not lead to confusion. 

To obtain ([50| . we, in fact, added to the right-hand side of ([TT]) the value 
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Let the interpolant p = p (x) be such that p^ = v". li x ~ 1 and = pj, then, 
by virtue of (|32|) -(|37 | as applied to the ceU [xi-i,Xi+i] x [i„,t„+i], wc find that 

E,: = -^Pr + 0((A-f)- (54) 

Thus, when we add of (|55)) to the right-hand side of (|47p . we, in fact, add the 
negative fourth-order diffusion term of (l54t to the right-hand side of the first 
differential approximation (1521) . Such a term, in contrast to a negative second- 
order diffusion term, stabilizes the system and produces dissipative effect (see, 

e.g., mi)- 

li x = 1 and p^ is estimated by the central-difference derivative, i.e. = 
(Pi+i - Pi-i)/(2Aa;), then 



{Axf 



6 

In such a case we obtain, instead of (1541). that 



p'r + 0{{Axf]. (55) 



\4 



E. = -^^Pr + 0((Ax)^j. (56) 

We can see from ([56|) that the estimation of the derivative p[ by the central- 
difference derivative leads to the same order of accuracy as it was in the previous 
example when di = p^. Furthermore, the scheme ()50p will be more dissipative 
in such a case, since the coefficient of the fourth-order diffusion term in ([55]) is 
three times bigger than the similar coefficient in (l54l) . 

If < >ir < 1 and p^ is estimated by the central-difference derivative, then 
we find 

E. - (1 - >c) ^pr + (1 - 4x) ^Pr + O [(Axf) . (57) 



Note that the estimations ([54|) . ([56|) . and (f57| will be valid if the derivative di 
need not be modified, i.e. (j42|) will be valid without implementing the algorithm 
Otherwise, there is a risk to add negative second-order diffusion terms 
to the right-hand side of (l52t . Let us consider, for instance, the case when 
K = 1 and the derivative p^^^ is estimated by the right- difference derivative, i.e. 
di+i = (pi+2 — Pi+i)/Ax, however p[_i is estimated by the central-difference 
derivative. In such a case we obtain: 

E. - -^-^P'U. - i^^P::% + O {(Axf) . (58) 

In view of ([551) and ([57)1 . we conclude that the part played by the parameter x 
in suppressing spurious oscillations produced by negative second-order diffusion 
terms can be very important. 

Interestingly, there is a possibility to improve the scheme (|45|) by introducing 
an additional positive numerical viscosity, i.e., using the vanishing viscosity 
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method [15], [32], such that the scheme's order of accuracy would increase up 
to 0((At)^ + (Ax)^). Such an approach [5] leads to the following second order 
central scheme 



^+0.5 



0-5 (.v,.+i + vj - K— (d^+i - d J - ^ 



f(v?) 



Ax 



8Ax 



(?u' 



(59) 



where df is the derivative of the interpolant at x = x^. The scheme (|59)) is 
abbreviated in [8] as C0S2. 

Let us develop a nonstaggered central scheme that will be of second order. 
Using the central differencing, we approximate (|22p at the point x = Xi, t = 
tn+05 by the following equation: 



df 



v"+i = v" - At ^ 



(9x 



Using Taylor series expansion, we approximate 



X=Xi, t = t„ + 0B 

n+0.5 



By virtue of the PDE system. 



9f(vr 



dt 
we find 



(60) 

in (|46l) with the accuracy 
(61) 



^ + O (At^) 



5f 9f 9u 

dt~ du' ~dt 



dt di dn 
dvL dvL dx 



dn 

dx 



(62) 



By virtue of (|49| and (16T|)-(|62|), we find 
vf_i +2vr+v7+i Ax 



ri+l _ "t-l 



9f 



-(dIV,-dILO-At^ 



a;— Xi, t—tn 



{Aty d_ (j^2^ 

2 9x V 9x 



5u 5u 



(63) 

Then, approximating the last two terms in the right-hand side of (|63p with the 
second order accuracy, we obtain 



'i-i 



2vr 



At 

2Ax 



(ff+i-fr-i) 



2 (Ax) 



2 [^"+0.5 



'i+1 



D 



i-0.5 



« - vr-i)] ' ^ = const > 0, (64) 



where T>^+Qr^ 



0.5 



(AD' 



(ADiJ 



the parameter is introduced by anal- 
ogy with K in (l45t . Scheme (|64l) coincides with the scheme MACl, (|50)) . pro- 
vided that <; = 0. The last term in right-hand side of ([M)) can be seen as the 
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non- negative numerical viscosity introduced into the first order scheme ([50]) . 
Thus, we are deaUng with the vanishing viscosity method [18], [32] and, hence, 
in view of [HI Theorem 3.3], the scheme, satisfies the entropy condition. 

To make this point a httle bit more clear, we consider a viscous perturbation of 
the system ([^^ . We associate with (1^^ the following parabolic system 



9Ue 



d_ 

dx 



f(Ue) 



d_ 

dx 



9uc 

dx 



e > 0, 



(65) 



where the right-hand side can be viewed as a viscosity term. It is assumed 
that Mock's assumption of admissibility [TS] p. 32] is valid, i.e. the matrix 
U" ■ is positive-definite. Here U = U (u^) denotes a strictly convex entropy, 
U" denotes the Hessian matrix of U. We recover solutions of (|22p as the limits 
of solutions to (|65l) as e — )• 0. If we divide the scheme COSl, (|50l) . by At and 
group all the terms of this scheme together in its left-hand side, then this group, 
in view of the above, is the finite difference approximation of the left-hand side 
of Eq. (p5|) on the cell [xi-i,Xi+i] x [tn,tn+i]- The right-hand side of Eq. (p5|) 
is approximated as following 



d_ 

dx 



dx 



D 



i-l-0.5 



-0.5 



(Ax) 



(66) 



where T>f±o,5 



0.5 



(Ar±i; 



If e = ^At/2, then we obtain the 

scheme (|64l) as the approximation of Eq. (j65|) on the cell [x^-i, x^+i] x [i„, i„+i]. 
Thus, owing to the last term in the right-hand side of (|64|) . the grid function 
v" in (|64|) can be considered as Lipschitz-continuous. Moreover, owing to this 
term, the scheme (|64)) is 0{{Ax) + (At) ) accurate provided that k ~ q — 1. 



The nonstaggered central scheme (|64|) . being a monotone approximation of the 
1-D equation ([22l) with the second order, will be abbreviated to MAC2. 



3.3 Stability of the developed schemes 

It is proven in ^ that the scheme C0S2, ([59]) . will be, in general, stable if 

u-^c2|H + a<i, a = ^^. (67) 

' ' Ax 

Let us note that the stability condition ()67|) is not exactly correct, namely, 
there exist infrequent situations when (|67p is valid nevertheless the scheme is 
not stable provided H = const. In particular, there could be a grid node 
(x,,i,*) where the following inequality must be, but not valid: 

<;HC2 + a<l, a-^^^. (68) 

Ax 

However, there is a possibility to take the parameters >c — >c{x,t), <j — (^(x^t), 
and H = ^{x,t) such that the scheme C0S2, ((SS]), will be stable under (|57)) . 
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For instance, if we take >c,(;,Cr = 1, H = 0.5 at ^ and H = at 

{x,t) = then ([68)1 will be also valid. 

In this subsection we will use analogous technic [5] to prove the stability of 
the developed nonstaggered central schemes. In view of Theorem[Hl the stability 
of the scheme MAC2, will be investigated on the basis of its variational 

scheme. It is assumed that the bounded operator A {— di (u) /du) in is 
Frechet-differentiable on the set J7u C M*^, and its derivative is bounded on Q,^. 
Considering that v" in (|59p is Lipschitz-continuous, we write 



v"||2 < CuAa:, Cy = const. 



'i+l ~ "i ||2 - ^u^-^; — ■-i^'i.cx,. (69) 

the second term in right-hand side of (j64p can be written in 



By virtue of 
the form 

(diVi - c-i) = f [Br • K+i - vr) - K-i ■ - v,Li)] 

Then, the variational scheme corresponding to ()64p is the following 
5v,r+i-0.25 (<5vr+i + 2<5vr + <5vr_i) 
At 



(70) 



2Aa 



(A^i . ^vIVi - Ati • <5v^i) 



+- [Br • (<5vr+i - <5vr) - a^.^ • (<5vr - s^u)] 



2Aa;2 
(Ai)= 



2Aa 



By virtue of (|42j) and ([231), we find 



< 4H, 



[(^DiVo.s) • (vr+i - vr) - (<5Dr_o.5) • (vr - vr_o] 



- [(<5Br) • (vr+i - vr) - (<5Ar_i) • (vr - vr_i)] 



(Ar)^ 



< A 



Thus, we may write that 

II^Dr+o.5LJnDr-o.5]L<2A: 



2 

max ■ 



(71) 



(72) 



(73) 



Then, by virtue of (12T|), (|69l), (172|)-(173]), and since < ^f, H < 1, we find the 
following estimate for the right-hand side of (|7ip: 



'2Ax2 



I (<5Dr+o.5) • (viVi - vr) - (<jDr_o.5) • (vr - v 



') • (vivi - vr) - {^-i) ■ (vr 



< 
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(ll'5Dr+o.5lU|vr-M - vrii, + ||<sD,_o.5ii, 

1 

8 



i {\\sMn2 ■ l|vr+i - vrii^ + \\sA-_,\\^ ■ ||vr - vr^jj < 



(l + 2C,2)C„MoAt, a = (74) 

Since the uniform stability with respect to the initial data implies the stability 
of scheme [JHl PP- 390-392] (see also Sec. 5]), we conclude, in view of (ffi)) . 
that the scheme ((7T|) will be stable if the following scheme (i.e. ((7T|) without 
the right-hand side) will be stable 

Sv^' = 0.25 (<5vIVi + 2Sv- + S^r^^) - ^ (A^+i • Jv^+i - A^^ ■ Sv^) 

- 1 [Br • (-jv^Vi - '^vr) - A^Li • (<5vr - ^v^lo] • (75) 

Since the bounded operator A (= di (u) /du) in ([22]) is Frechet-differentiable 
on the set Vl^ C M*^, and its derivative is bounded on Q,^, the operator A is 
Lipschitz. In view of (|69|) . we can write that || (At/Ax) (Af_j — A-':^^) • 5vf || < 
C,At\\5^^^^\\, Ca = const. Thus, adding the term (At/Aa;) (Ar_i - Af+i) • 
(5v" to the right-hand side of ([75]) has no effect on the stability of ([75]) . To 
investigate the stability of (TfSt we rewrite it in the following form, where the 
term 0.5(Ai/Ax) (Af_i - A"}^^) ■ Jv^ is added. 

+ (0.51 + E^,_i + E^,;+i) . + (0.251 - E^.+i) • Jv^Vi, (76) 

where 

= - + (78) 

We write, in view of ([23|) and (|42|), that the spectrum s (E-^) C [£min,i, ^max.i] ' 
where j — i±l and 

^min.. > -0.5?C2 - 0.5a, imax,^ < 0-5^^ " 0.25^^^ + O.SC^- (79) 

Hence, by virtue of |71 Theorem 2.10] we find that the scheme ((76)) will be stable 
if 

max (10.25 -AK 10.25 A|) 0.5, Vz,n. (80) 

^^[-C'min.i'-E'max.i] 
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We obtain from (150)) the following sufficient conditions for the stability of the 
variational scheme (iTlT) 

- 0.5<;C2 - O.Sa > -0.25, (81) 

0.5xH - 0.25<^C2 + 0.5a < 0.25, C,. = (82) 

Ax 

Let us note that the value of ij^^^x i G3) can rarely be attained. It could be the 
case at the point (xi,t„) only if ||B"||2 = 4H, ||A"^^||2 — Amax, and ||A"||2 — 
(or, what is the same, jjA^.J^ = 4H, ||A"_;^||2 = Amax, and |iAf||2 = 0). 
Moreover, considering the modification of the scheme (|76l) . where D"_q g is 
replaced by (A"_]^)^, namely 

E".-. = fA:'.-.^(A?_,)'-|^Ar_,. (83) 



ES.« = |B?-<i^(A?«)^ + 5|:A?«, (84) 



we note (see Proof in 7; Theorem 2.10]) that the stability of the modified 
scheme, jTSD, ([Ml) -(El), implies the stability of the scheme d75)) - ((75)l . since the 
operator A is Lipschitz. Hence, instead of ([5^ . we will consider less rigid 
requirement: 

0.5>fH - 0.5^C2 + 0.5C,. < 0.25. (85) 

Let = 1, then, by virtue of ([5T|) and (|55|) . we find the stability condition for 
the variational scheme ([7T|) : 

Cr < 0.5 (\/3 - l) , >^H < 2 - V3. (86) 

Thus, in view of Theorem H the scheme MAC2, will be stable if §^ will 
be valid. 

Notice, if <j -> in ([5T|) and ([5^ , then we obtain the stability condition for 
the scheme MACl, ([501 : 

>^H + a<0.5, a = ^^^. (87) 

It is significant that the stability conditions (ISTI) . (|55|) . and hence (|5S|) . is 
found by virtue of |7l Theorem 2.10] and, therefore, it is assumed that Ef^ is 
Lipschitz, i.e. E"j_|_^ = E"j_^ + 0(Aa;). For the sake of brevity, we associate the 
Landau symbol 0{Ax) with Lipschitz-continuity of a grid function whose incre- 
ment can be estimated in norm by a constant times Ax for Aa; small enough. 
Note that E"^ is Lipschitz if A", A", and B" are Lipschitz. It is easy to see that 
A" will be Lipschitz if u" (x) (the solution to (|65|) ai t = tn) will be Lipschitz- 
continuous. Actually, since the bounded operator A (u^) (= (9f/9u)|^^^ ) is 
Frechet-difFerentiable on the set il^ C M^^, and its derivative is bounded on 
fiu, we write: ||Af_^i "^"11 ^ ^"^ ~ - CaCuAx, where Ca, C„ 
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= const. We may guarantee that A" and B" will be Lipschitz if u" (x) will be 
differentiable. Actually, we may write, by virtue of p8l) . that = I + 0{Ax), 
and Bi = I + 0{Ax). Hence, Lipschitz-continuity of A" and Bf is obvious. 
However, if u" (a;) is only Lipschitz-continuous, then we can not guarantee that 
A" and B" are also Lipschitz-continuous. In such a case a generalization of 
[3 Theorem 2.10] can be useful, namely the case when the scheme coefficients 
depend on several matrices. Then we obtain that (|76p . ([83|) - (l84|) will be stable 
if 

max (|0.25-Ai| + IO.25-A2I + |0.5+Ai + A2I) s$ 1, Vi,n. (88) 

By virtue of we find that the scheme MAC2, will be stable under the 
same conditions ([5T|) . (I55|) . and, hence, under (1551) . 

Let us note that the stability conditions, are, in general, apt to be 

unduly rigid. Assuming that H = H(a;i,i„) as well as Cr = Cr{tn), and using 
Weyl's inequalities (see, e.g., [22l Chap. 3]), we can find less rigid stability 
conditions than Let A"'*^ denote the k — th singular value of A", A^^^^^ ^ = 
maxA"''^, Ajjja^x — '^^^'^max i; ^'^'^ o^'^ dcnotc the k — th singular value of Af , 

k i ' 

•^min i — niina"''^. Then, given the pre-assigned H = H = const and Cr = Cr = 

const < 0.5, we find the sought-after stability conditions at i = t„, employing 
the following two-step procedure. 

Algorithm 10 

Step 1. Given Cr ^ Cr < 0.5 and t" = AxCr/X^i,^, we, by virtue of Weyl's 
incquahties [331 Chap. 3], estimate U^^^f 

i".. , < 0.5>^H" - 0.5 ( (A"., y + 0.5^A" . 



Ax J AX 



In view of dHS]) and §9\j, we find 
H" = min 



W f_l (^max,i) t" 1 

' ' Ax / X xAx 2>^ 



(90) 



Step 2. Given Hf, (|90p . wc find the derivative d" using (l43l) and, by virtue 
of (|35|) . the diagonal matrix A", and hence a^^^^. We, by virtue of Weyl's 
inequalities [22, Chap. 3], estimate m^^^ .f. 

iJJ.i„,.>^^-0.5(^— j (a;;,,,)'~0.5^A;;,.,,. (91) 

Let r" be used instead of r" in ([9T|) . In view of (|88|) and (|9T|) . the following 
inequality must be valid. 



.. , 2 

r 



2 t" 

A ; (^max,i) ~ X~'^max,i ^ -0.5. (92) 

Ax / Ax 
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We obtain from (1^^ : 



Aa 



max.t ^ y ' inin,t (93) 



Thus, given r" (= Aa;C'r/A^^3^^) and r", (|93|) . we find the time increment, 
Ai", at t = i„: 

At" = min (T",minT^j ^ C;^ = ^ J^^" . (94) 

4 Operator splitting schemes 

By virtue of the operator-splitting idea [13], [T^, [32] (see also LOS in [H]), 
the following chain of equations corresponds to the problem (|T]) 

^^ = -q(U), t„<i<t„+o.5, U(x,t„) =U"(x), (95) 
2 at T 

1 ^TT ^ a 

a^fj- (U) = 0, < t < i„+i, U (x,t„+o.5) = U"+"-' (x) , (96) 

where U" (x) denotes the solution to (IMl) at i = t„, U"+°-^ (x) denotes the 
solution to (|95l) at i = t„+o.5. If a high-resolution method is used directly for 
the homogeneous conservation law (j96p . then it is natural to use a high-order 
scheme for (IMt . As applied to, in general, stiff (r ^ 1) System (jlj, the second 
order schemes can be constructed on the basis of operator-splitting techniques 
with ease if ()95p will be approximated by an implicit scheme and (|96p by an 
explicit one, see Proposition 4.2 in 7 . As an example, let us develop a central 
scheme for a 1-D version of (jlJ. After operator-splitting, the 1-D equation can 
be represented by the chain of equations, namely (|95l) and 

1 (9U 8 

2 ^ + (U) = 0, t„+o.5 < t < t„+i, U {x, i„+o.5) = U"+°-5 {x) . (97) 

Let us first consider the case when the following first-order implicit scheme be 
used for ^ 

^-q(vr-^), (98) 



and a central scheme with nonstaggered grid cells will be used for ([^7]). We 
rewrite the scheme MACl, ((5(I]) . to read 

vr+i = 0.25 (vl!^-^ -f 2vr°-^ + v^!+"-^) - (d^!+"-' - dl!^-') 

At .^„+^o.5 _ ^ ^ ^ ^^^^^^^ < >^ < 1, (99) 



2Ax 
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where d""*""'^ denotes the derivative of interpolant at a; = Xi. It is clear that the 
scheme ([M)) approximates (P7|) with the accuracy 0{At + (Ax)^), however, in 
view of Proposition 4.2 in \7], the scheme (IM|) - (IM1) . taken as a whole, is of the 
second order approximation for the 1-D version of ([T]). 

Let us develop another nonstaggered central scheme approximating a 1-D 
version of (H]) with the accuracy 0((Ai)^ + (Ax)^) and such that its compo- 
nents (after operator splitting) will be of the second order. It can be done on 
the basis of the second order scheme (|99p with ease. Notice, adding to 

and subtracting from Equation (4.10) in 7_ (rewritten for tn < t < i„+i) the 
same quantity is equivalent to adding this quantity to (|99)) and subtracting it 
from dMl). Let 0.125 (At)^ {d'^\]/dt^)'l^^-^ be this quantity, then we obtain the 
following scheme, instead of (|98l) . (l99l) . 



n+0.5 ^ 



n+0.5 



n+0.5 



(100) 



+1 = 0.25 (vf+f^ + 2v: 



n+0.5 I , n+0.5\ (-jri+O.S 



2Aa 



m+0.5 



A 

(Ai)^ ^ a^u 



(101) 



Thus, the scheme ()100p as well as the scheme (llOip are of the second order, and 
the scheme (|100p - (jl01l) . taken as a whole, is of the second order as well. Using 
Taylor series expansion, and central differencing, we find 



1 /At 

2 It 



A \ dt 

n+0.5 



n+0.5 



At fdV 



2 \ dt 

We obtain, by virtue of ([62]), (HZl), that 



ri+0.25 



0{{Atf), 

o{{Aty 



(102) 
(103) 



dHJ_ _ d_ /dr\ _ d_ /di\ _ d_ / 2 9U 
'W ^ ^ di [d^ J ~ [dt J fel 'a^ 



(104) 



where A =9f /9U. Then 



" d 




n+0.5 ^ 








dx 




" A^ 



pin+0.5 ^i+l ^' 



n+0.5 



Ax 



,n+0.5 



7~>n+0.5 
^1-0.5 



n+0.5 



Aa; 



O (Ax)' 



(105) 
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where B^^+^f = 0.5 ((Afi" ')' + (Af+°-5)2). By virtue of m, we 
rewrite Scheme (|100l) - (|10ip to read 

n+0.5 _ ,,n I ^^„n+0.25 



V, = V, 



— qr+"-'', (107) 



T 



c(At) r_„+o.5 / n+0.5 ,,n+0.5\ i~,n+0.5 (n+0.5 _ , n+0.5\] 

2 (Ax) 



2Aa; 



If = 0, then ([T08| coincides with (HH), being 0(At + (Ax)^) accurate, li k = I 
and = 1, then Scheme (|106p - (ll08p approximates the 1-D version of (H]) with 
the accuracy 0((At)^ + (Aa;)^). The scheme (jlOSp coincides, in fact, with the 
scheme MAC2, dM]), and hence the scheme (jlOSp will be stable if dHU) will be 
valid. Notice, the less rigid conditions for the stability of the scheme (|108p can 
be found by means of Algorithm [101 

Let us note that in practice (e.g., [32], [48]) the operator-splitting techniques 
find a wide range of application in designing economical schemes for Eq. (|96p 
in the domains of complicated geometry. The resulting method, in general, will 
be only first-order accurate in time because of the splitting [32], [48]. Thus, in 
line with established practice we will replace the multidimensional Eq. (|96p by 
the chain of the one-dimensional equations: 

2N~dt' ^ dx]^' ^'^^'^ ^' * - Wo.5+i/(2JV), (109) 

where Uj {x,tn+0.5+{j-l)/{2N)) = '^j-l (2^, ^n+0.5+(i-l)/(2A')) , j = 1,2, . . .iV, 

Uo (x, in+o.s) denotes the solution to (|95l) a.t t = tn+o.5- Eq. ([95]) will be 
approximated by the first-order implicit scheme, (|98p , or a second-order implicit 
Runge-Kutta scheme. In particular, it can be used (|106p - (jl07l) or the following 
implicit Runge-Kutta scheme 

n-l-0.25 _ n-l-0.5 _ ^„n+0.5 n+0.5 _ , ^ n+0.25 (^^n'\ 

2t ' i i ' J. ' V-*^-*^^/ 

since these schemes possess a discrete analogy to the continuous asymptotic 
limit (see, e.g., [53], [37], [H])- Let us note that the scheme (|110p is an implicit 
analogue to the well known Runge-Kutta scheme, which is referred, originally 
due to Runge, as modified Euler method [5D]: 
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5 Examples 



In this section, we are mainly concerned with verification of the second order 
central schemes C0S2, (|59|) and MAC2, (l64t . as well as the splitting schemes 
p^ - pUg]) . and pTO)) . 

5.1 Scalar non-linear equation 

As the first stage in the verification, we will focus on the following scalar 1-D 
version of the problem ([T]) : 

du d 

^ + ^/H-o, -T e M, < t < r,,ax; u{x,t)\^^^ = u"{x). (ii2) 

Let us first compare the schemes MAC 2, (|64l) . and C0S2, (f59|) . and demon- 
strate that these schemes are of the second order. We consider the linear 
transport equation, i.e. (jll2l) with / (u) = u, subject to the initial data: 
M° (a;) = sin(7ra;). The numerical solutions were computed under k = <, — 1. 
The scheme C0S2, will be stable, in view of (jSZl), (jMl), if we take the CFL 
number Cr = VS- 1, and H = 2 - VS. The scheme MAC2, ([Ml), will be stable, 
m view of dMl), under Cr = 0.5(V3 - 1) and H = 2 - V3. We will also verify the 
scheme COSl, (gS]), under Cr = 0.5 and H = 0.25, as well as the scheme MACl, 
((50l). under Cr — 0.25 and H = 0.25. Notice, in the case of schemes C0S2 and 
COSl the CFL numbers are two times higher than in the case of schemes MAC2 
and MACl, respectively. The reason is that the staggered schemes, C0S2 and 
COSl, are solved twice during the time increment. At, in contrast to the non- 
staggered schemes MAC2 and MACl. Li errors, at t = 10, < a; < 2, versus 
the number of nodes are depicted in Table [1] 



Table 1: Li errors versus the number of nodes (N) 



N 


1280 


640 


320 


160 


80 


COSl 


2.4 X 10-^ 


4.9 X 10-^ 


9.8 X IQ-^ 


2.0 X 10-^ 


3.9 X 10-^ 


C0S2 


1.5 X 10-^ 


6.2 X 10-^ 


2.5 X IQ-* 


1.0 X 10-^ 


3.9 X 10-^ 


MACl 


2.5 X 10-^ 


4.9 X 10-^ 


1.0 X IQ-^ 


2.1 X 10-^ 


4.3 X 10-^ 


MAC2 


5.6 X 10-^ 


2.4 X 10-* 


9.8 X IQ-* 


3.9 X 10-^ 


1.5 X 10-^ 



Note (Table [Ij that Li errors in the case of the scheme C0S2 are approxi- 
mately coincide with those appeared in the case of the scheme MAC2 under the 
double number of nodes. It is associated with the fact that in the case of the 
staggered scheme C0S2, (|59|) . the space increment is, actually, two times less 
than in the case of the nonstaggered scheme MAC2, (|64p . 

Now we will solve the inviscid Burgers equation (i.e. / (u) = u^/2) with 
the following initial condition 

«(x,0) = I ""J"' llll'^'l''] , hR>hL, uo^ const ^0. (113) 
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The exact solution to (|112p . (|113p is given by 



where T — 2S / uq, S = Hr — h^, 

\ f^"o, Hl <x<b, b^uot + hL 
ui{x,t) = < Uo, b < X < O.duot + Hr , (115) 

[ 0, X < Hl or X > O.buot + 

{2S{x-hL) I ^ ^ T 

U, X < tiL or X > L 

L = 2 + Q.buoS {t-T) + hL. (117) 

First, it wih be used the first order in time schemes MACl, ([5U|. and COSl, 
(|45p . The numerical solutions were computed on a uniform grid with spatial 
increments of Aa; = 0.01, the velocity Mq = 1 in (jll3l) . ft-L = 0.2, Hr = 1. In 
view of the stability condition, (I57|) . for the scheme MACl, we take the CFL 
number Cr = 0.25, the monotonicity parameter H = 0.25, and the parameter 
K = 1. Having regard to the stability condition (j67p under = 0, we verify 
the scheme COSl, (gS]), under Cr = 0.5, H = 0.25, and k = I. The resuhs of 
simulations are depicted with the exact solution in Figure [21 




Figure 2: Inviscid Burgers equation. The schemes COSl (a) at given Cr = 0.5 
and MACl (b) at given Cr = 0.25 versus the analytical solution. Crosses: 
numerical solution; Solid line: analytical solution and initial data. Aa; = 0.01, 
H = 0.25. 

We note (FigurelJ) that the schemes COSl, (gS]), and MACl, O, exhibit 
similar results in spite of the fact that in the case of the scheme COSl the space 
increments are, virtually, two times less than those in the case of MACl. It 
should also be mentioned the absence of spurious oscillations in the numerical 
solutions (see Figure [2]). However, if we take H = 0.5, then the spurious oscilla- 
tions can be produced by the scheme COSl (see [HI P- 2804]). Using Algorithm 
[TOl at given H = 0.5, in the case of scheme MACl we obtain very slight but 
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significant spurious oscillations, which indicate that the boundary maximum 
principle is violated. For instance, the maximum positive value of the depen- 
dent variable, w, at i = 1 is 0.4% above the maximum value oiv at the boundary 
t = 0. These spurious oscillations can be eradicated by decreasing the parame- 
ter >c. Particularly, the spurious oscillations disappear ii >c = 0.7, however, this 
introduces an additional numerical smearing. The results of simulations are not 
depicted here. 

To test the schemes MAC2, dH]), and C0S2, the inviscid Burgers 

equation was solved under the initial condition (|113p . The numerical solutions 
were computed under the same values of parameters as in the case of the schemes 
MACl and COSl, but Cr and H. In view of the stability condition, ([55]) . for 
the scheme MAC2, we take the CFL number Cr ~ 0.5(-\/3— 1). It is also used 
Algorithm [TU] at given H = 0.5. The scheme C0S2 will be stable, in view of 
(p7|) . if we take Cr — \/3 — 1 and H = 0.5. The results of simulation are 

depicted with the exact solution in Figure [31 




Figure 3: Inviscid Burgers equation. The schemes C0S2 (a) at g 
73 -1 and MAC2 (b) at given Cr = 0.5 (n/3- 1) versus the analytical solution. 
Crosses: numerical solution; Solid line: analytical solution and initial data. 
Ax = 0.01, H = H = 0.5. 

We note (Figure [3]) that the schemes C0S2 and MAC2 exhibit a typical 
second-order nature without any spurious oscillations. The scheme MAC2 ex- 
hibits more numerical smearing than the scheme C0S2. Such a phenomenon 
has, at least, two reasons. First, in the case of the staggered scheme C0S2 
the space increment is, in fact, two times less than in the case of the nonstag- 
gered scheme MAC2. Second place, using Algorithm [TU] at given H = 0.5 we, 
in general, obtain H" < H at the grid node {xi,tn), resulting in an additional 
numerical smearing. 



5.2 Hyperbolic conservation laws with relaxation 

Let us consider the model system of hyperbolic conservation laws with relaxation 
developed in [45] : 
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dz d 1 

- + -az = -Q(u;,z), (119) 

where 

Qiw, z) = z — m{u — uq), u = w — qqz, (120) 

T denotes the relaxation time of the system, qo, m, a, and uq are constants. The 
Jacobian, A, can be written in the form 

^ 1^ w - qoz + a -go (w - qoz) | ^^^^^ 

The system (|118p - (jll9p has the fohowing frozen '45' characteristic speeds Ai = 
a, X2 = u + a. The equilibrium equation for (|118p - (|119p is 

where 

TTl 

= w — qoz^,, = {w — Uq) . (123) 

1 + mqo 

The equilibrium characteristic speed A* can be written in the form 

A, (w) - Jfll!^ + a. (124) 
1 + mqo 

Pember's rarefaction test problem is to find the solution {w, z} to (IllSp - 
(|119|) . and hence the function u — u{x,t), under r — >■ 0, and where 

{w,z} = l {^L,z*in^L)h x<XQ 

< UL ^ WL - qoZ* (wl) <ur^wr- qoz^ (wr) . (126) 

The analytical solution of this problem can be found in [45 . The parameters 
of the model system are assumed as follows: go — ~1, ™ — ~1, uq = 3, a = 1, 
T = 10~*. The initial conditions of the rarefaction problem are defined by 

ul = 2, =^ zl = m {ul - Uq) = 1, wl = ul + qaZL = 1, (127) 

Ur = 3, =^ Zr = m {ur - Uq) = 0, Wr = Ur + q^ZR = 3. (128) 

The position of the initial discontinuity, xq, is set according to the value of a 
so that the solutions of all the rarefaction problems are identical [45]. Let a 
position, x]^, of leading edge or a position, x^, of trailing edge of the rarefaction 
be known (e.g., = 0.85, a;^ = 0.7 in IS]), then 

^0=4- (tT^ + a)t^xi- + a) t. (129) 

V 1 + mqo ) V 1 + "^90 / 
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At t = 0.3, under (IT?f|) -(in5 1) we have gS] 



X < 0.7 

0.7 < X < 0.85 . (130) 
X > 0.85 

The results of sumilations, based upon the schemes MAC2, (|64l) . and C0S2, 
([551) . together with (jllOp . under different values of a grid spacing {Ax = 10~^, 
5 X 10^**) are depicted in Figure S] In view of the stability condition, ((86| . for 
the scheme MAC2, we take the CFL number Cr = 0.5(a/3— !)• It is also used 
Algorithm [TU] at given H = 0.5. The scheme C0S2 will be stable, in view of 
(P7|). (PI) , if we take Cr = 73 - 1 and H = 0.5. 



2 + 



2;-0.7 



0.85- 

3, 



0.7' 




Figure 4: Pember's rarefaction test problem. The schemes C0S2 (al, a2) at 
given Cr = VS- I and MAC2 (bl, b2) at given Cr = 0.5 (%/3 - l) versus the 
analytical solution for u. Dashed line: numerical solution; Solid line: analytical 
solution. Time t = 0.3, H = H = 0.5, (al, bl): Ax = 0.001, (a2, b2): Ax = 
0.0005. 

One can clearly see (Figure 2]) that the schemes MAC2 together with the 
Runge-Kutta scheme (jllOp as well as C0S2 together with (IllOp are free from 
spurious oscillations. Recall that in the case of the staggered scheme C0S2 the 
space increment is, in fact, two times less than in the case of the nonstaggered 
scheme MAC2. Moreover^ using Algorithm [TOl in the case of scheme MAC2, we, 
in general, obtain H" < H at the grid node {xi,tn), resulting in an additional 
numerical smearing. Nevertheless, the scheme MAC2 can exhibit (Figure H]) 
even less numerical viscosity than the scheme C0S2. 
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5.3 1-D Euler equation of gas dynamics 

In this subsection we apply the second order schemes C0S2, ([5^ . and MAC2, 
([M]). to the Euler equations of gamma-law gas: 

^^^S^ + (u) = 0, xeR,t>0; u{x,Q)^u"{x), (131) 
ot ox 

u={ui,U2,U3} ^{p,pv,e} , F (u) = {pv,pv +p,{e + p)v} , (132) 

e = — 1 — pv'^, 7 = const, (133) 

7 — 1 2 

where p, v, p, e denote the density, velocity, pressure, and total energy respec- 
tively. We consider the Riemann problem subject to Riemann initial data 

"°(-)-{u^ ItZ ^ -^^-n = const. (134) 

The analytic solution to the Riemann problem can be found in [321 Sec. 14]. 

We solve the shock tube problem (see, e.g., [S], [53], [52], [33]) with Sod's 
initial data: 









r 0.125 ] 


Ul = 1 














[ 0.25 J 



(135) 

Following Balaguer and Conde 5 as well as Liu and Tadmor |33j we assume 
that the computational domain is < < 1; the point xq is located at the 
middle of the interval [0, 1], i.e. xq — 0.5; the equations (|13ip are integrated up 
to t = 0.16 on a spatial grid with 200 nodes as in [5] and in [3^. Aiming to 
compare the schemes C0S2 and MAC2, we, in view of the stability conditions 
([67|). ([68|. and take the CFL number Cr = ^3 - 1 (the scheme C0S2), 

and Cr = 0.5(%/3- 1) (the scheme MAC2). In both cases we take >fH = 2 - \/3, 
where >r = 0.94. The reason why the parameter >c has to be less than unity 
is discussed in Section [3?^ (see also |8] Sec. 4]). The results of simulations are 
depicted in Figure [5](left column) and Figure [SK left column). 

The results using the scheme C0S2 (Figure [S] left column) are not worse 
in comparison to the corresponding third-order central results of [33l p. 418] 
as well as to the results obtained by the fourth-order non-oscillatory scheme 
in [5] p. 472]. The simulations in [5] and [S^ were done under At = O.lAx 
(i.e. 0.13 < Cr < 0.22). The fourth-order scheme [SI p. 472] gives a better 
resolution but, in contrast to the scheme C0S2 (Figure [5l left column), can 
produce spurious oscillations. Let us note that the second-order scheme C0S2 
can give results analogous to those obtained by the fourth-order scheme [S] p. 
472]. Namely, the higher will be the values of >€ and H, the better will be 
the resolution (see, e.g.. Figure [5l right column). However, if >f = 1, then the 
scheme C0S2 will be oscillatory [H p. 2816]. It is apparent that x should be 
reasonably less than unity to suppress spurious oscillations, while it should be 
as close to unity as possible to lose less in the order of accuracy of the scheme [8] 
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Figure 5: Sod's problem. Time t = 0.16. The scheme C0S2 versus the analytical 

solution. = 2 - yc = 0.94 (left column) and H = 0.35, k ^ 1 (right 
column). Crosses: numerical solution; Solid line: analytical solution. Cr = 
V3-1, Aa; = 0.005 
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Figure 6: Sod's problem. The scheme MAC2 versus the analytical solution. 
Time t = 0.16. Cr = 0.5(\/3 - 1), ^ 2 - y^, x = 0.94. Crosses: numerical 
solution; Solid line: analytical solution. Ax = 0.005 (left column), Ax = 0.0025 
(right column). 
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p. 2816]. The results of using the scheme MAC2 (Figure [SI left column) are not 
worse in comparison to the corresponding second-order central results of [33l p. 
418], in spite of the fact in the case of staggered scheme the space increment 
is, in effect, two times less than in the case of nonstaggered scheme. For the 
sake of illustration, the results using the scheme MAC2 under Aa; = 0.0025 are 
depicted in Figure [S] (right column). It is easily seen that the scheme MAC2 
gives a better resolution than the scheme C0S2 (Figure [5l left column). 



5.4 3-D axial symmetric gas dynamics 



We consider an adiabatic expansion of a vapor-plasma plume, produced by 
a laser beam on a target plane surface, into vacuum [3], i.e., the so called 
Anisimov's problem. It is assumed [3 that the expansion of the plasma plume 
may be described by the ideal gas dynamic equations with a constant adiabatic 
exponent 7. Let the target surface be perpendicular to the axis z and located 
at z = 0. Taking into account the symmetry of the plume with respect to the 
axis z, the gas-dynamic equations can be written as follows. 



it ^^^^ 

Ft 

dpE 

pE = 



dp 1 d jrpvr) 
dt r 

1 d 



r dr 
d_ 

dz 



ld_ 

r dr 
P 

7-1 



dr 

rp (Vr) 

rvr {pE - 
+ O.bpv^ 



P) 



_^ d{pv^) 
dz 

oz 
1 d 

-— [rpv, 
r or 

dz 



-0, 



2 2,2 

V = v„ + v^, 



dp 
dr 
dp 
dz 



7 = const 





(136) 


0, 


(137) 


0, 


(138) 


0. 


(139) 




(140) 



Let i?o and Zq will be the initial radius of the plume at z = and the initial 
height of the plume in z-direction, respectively. It is also used in (3l the mass 
Mp of the plume and its initial energy Ep as the input data. We will use the 
following values as the reference quantities: = i?o, = ^ (57 — 3) Ep /Mp, 
= p^ = Mp/ [RqZo), Pt, = p^v^. It is assumed that the system 

(|136p - (|139p is written in the non-dimensional variables, since the dimensionless 
form of the gas dynamic equations coincides with (I136p - (|139p . up to notations. 
The pseudo- analytic solution to (|136p - (|139p can be written in the following form 



p{r,z,t) 



(7) 



9 2 

^2 ^2 



1/(7-1) 



, ^^Ht), v^vit), (141) 



p{r,z,t) = 



I2' (7) 



57-3 
Vr {r,z,t) 



j,2 ^2 



7/(7-1) 



Zo 

Rq 



r-, (r,z,t) 



^~ dt' dt' 



(142) 
(143) 
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where 



f2 



C (0) = 1, 77 (0) = a, U0) = V (0) = 0, (144) 



h (7) = 
h (7) 



TT 



3/2 



2 

3/2 



7 



2(7-1) V7-I 



7-1 
7 



7-1 



(145) 



(146) 



r(-) is the Gamma-function. Out of the plume, i.e. provided (r/^)^ + (z/77)^ 
> 1, the values of density {p) and pressure {p) as well as the components of 
velocity {vr,Vz) are equal to zero. 

Thus, the problem, (I136p - (|139p . is reduced to the system of ordinary differen- 
tial equations (ODEs) (|144p . The system of ODEs, (I144p . is solved numerically 
by the Runge-Kutta method with the adaptive step-size control [JS]. To check 
the accuracy of the calculations, the integral of energy [3] is used: 



^ +0.577^ 



7-1 



t2 



7-1 



Aiming to use the schemes C0S2, ((591), and MAC2, 
system p36l) - p39)) . we rewrite p^ - ((T39)) to read 



d_ 

dt 



dp 

dt ' 

d_ 

dr 
d_ 

dz 



d{pvr) , d{pvz) 



dr 



p{VrY 



p{Vz) 



dz 



pVr 

r 



P 



P 



^ ^ {pVzVr) = 

oz 



p{Vr 



r 

PVzVr 



dpE d_ 
dt dr 



[Vr{pE+p)] 



d_ 

dz 



dr 

[v, (pE + p)] 



Vr {pE + p) 



(147) 

for solving the 

(148) 

(149) 
(150) 
(151) 



The initial conditions are the following (in details, see (|141|) - (|146|) under 
t = 0): 



Vr =Vz=0, p = I 



1 / -> o\ 1/(7-1) 

(7) f 1 - - (z/cr) j , p/p^ = const. (152) 



The boundary conditions are established on the basis of the reflection concept 
[2]. At r = we assume that the axis z is a reflection line. It prohibits any 
normal flux of mass through the boundary r — 0, i.e. Vr = 0. Moreover, it is 
assumed that the pressure (p), density (p), and tangential velocity (vz) are even 
functions of normal distance to the axis z while the normal velocity {vr) is an 
odd function of r. It is also assumed that the plane z = is a reflection surface, 
i.e. the pressure (p), density (p), and tangential velocity (vr) are even functions 
of normal distance above the target surface while the normal velocity (vz) is an 
odd function of z. 
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By virtue of the operator-splitting method [5], [T3], [H], [22] (see also LOS 
in [48]), the system (|148p - (ll5ip may be approximated by the following chain of 
equations: 



Idu d ~ , . ^ , 



Idu 
2~dt 



-q(u), t 



n+0.5 



< t < t„+i, uL 



^ri+0.25 



-,n+0.5 



(153) 

(154) 
(155) 



where u ^ {p ,pvr, pi! z, pE} ,i{u) = {pVz,pVzVr,p{vzf+P,Vz{pE+p)} , 

U = {p , pVr, pVz, pE] , f (u) = {pVr , P {Vrf + P, PVzVr, Vr {pE + p)] , 

u = {p , pwr, pw^, P-Ej , q(u) = {pUr , p(ur) , pw^-Or, {;r(pi;+j5)| . 

We will use the schemes C0S2, and MAC2, fo solving Eqs. (ITCHj) - 
p54p . To solve Eq. (|155p it will be used the Runge-Kutta scheme (modified 
Euler method), (|llip . Let us note that every point on the axis r = is a singular 
point for Eq. (|155p . Assuming that all terms at left-hand side of Equation (I148P 
are bounded values at a vicinity of r = 0, we find that — > as r ^ and, 
hence, q(u)|^^Q ~ 0. Thus 



j.^q(u)U 

r— >0 r 



lim 



r>0 



q(fi)l.=0 5q(u) 



9r 



r=0 



Thus, we have at r = 0: 

1 9u 9q (u) 



2 9i 



dr 



U 



n+0.5 



(156) 



(157) 



Let Ar, Az denote the spatial increments and let = iAr, = jAz, 
i, j = 0, 1, 2, . . . . The Runge-Kutta scheme, for Eq. (|155p can be written in 
the form: 

-n+0.5 



-ri+0.75 
i-J 



At 



2r, 



vn-l-1 _ n+0.5 



At ,.„^ 

' 7, 



n+0.75\ 



z-1,2,..., j=0,l,2,... , (158) 

T 



Where v^^f = {(p)::f , {pVrY:f , {Pv.ty, ^PETIT ] ' ^ = O'^' 0-75, 1- If 
i = 0, then, in view of (|157p . the Runge-Kutta scheme is the following 



At aq 



--n+0.75 _ ,->Ti+0.5 



n+05 



r=0 



_ -n+0.5 



At^ 

or 



n+0.75 



, J > . (159) 



r=0 



Notice, we are using, for the sake of convenience, the same notation for the 
components of, in general, different vectors u and v. To derive the scheme at 
i = 0, we will use the reflection concept f2|. We introduce a node i = —1, i.e.. 
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r_i = — Ar. At this node, in view of the reflection concept, we have: v"]'j = 

Then, using central differencing, we obtain, by virtue of (I159p . the following 
difference scheme for Eq. (I157P : 

^ = v^r - ^ [q (-rr ) - q (v-tr )] ' 

= nr-' - ^ Kr-") - q (-"1°-")] , J = 0, 1, 2, . . . . (160) 

The equations p36p - (|139p are integrated up to i = 1 with a = ZqXRq = 
0.1. It is assumed for the schemes C0S2, and MAC2, ((Ml), that the 

monotonicity parameter H — 0.2 and the parameter >c ~ 1. In view of the 
stability conditions, (|67| we take Cr — 0.95 for the scheme C0S2. Following to 
[23], we take Cr = 0.475 for the scheme MAC2, in contrast to the sufffcient (but 
not necessary) stability conditions ([86]). It is assumed, in the case of scheme 
C0S2, that the spatial increments are the following: Ar = 10^"^, Az = 5 x 10"'' 
if < i < 0.1; Ar = 2 X lO'^, Az = 10"^ if 0.1 < t < 0.4; Ar = 2.0 x lO^^^ 
Az — 2.0 X 10^^ if 0.4 < t < 1. Taking into account the fact that in the case of 
the staggered scheme C0S2 the space increment is, actually, two times less than 
in the case of the nonstaggered scheme MAC2, we assume for the scheme MAC2: 
Ar = 5 X 10-'', Az 2.5 x 10"* ff < i < 0.1; Ar = 10-^, Az = 5 x 10"'' if 
0.1 <t< 0.4; Ar = 10"^, Az = 10"^ ff 0.4 < i < 1. The results of simulations 
as well as the analytical solution are depicted in Figures [3 [SJ El [TOl 

We observe (Figures [71 [3 [9l [10]) that the numerical results obtained via 
the schemes C0S2 and MAC2 are practically coincide. We also observe that 
the relatively large deviations between numerical and analytical solutions are 
occurred in the area of the relatively large second derivatives (i.e. their absolute 
values) of the primitive variables over the spatial coordinates. Hence, these 
relatively large deviations may be occurred in the area of small values of the 
variables in the vicinity of the front (see Figures [7] [S] [9] [T0| as well as in the 
area of relatively large values as well (see Figure [8]). 
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0.95) and MAC2 (Cr = 0.475) schemes versus tlie analytical solution, a = 
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